{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def load_data(filename):\n",
    "    \"\"\"载入数据。\"\"\"\n",
    "    xys = []\n",
    "    with open(filename, 'r') as f:\n",
    "        for line in f:\n",
    "            xys.append(map(float, line.strip().split()))\n",
    "        xs, ys = zip(*xys)\n",
    "        return np.asarray(xs), np.asarray(ys)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 不同的基函数 (basis function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def identity_basis(x):\n",
    "    ret = np.expand_dims(x, axis=1)\n",
    "    return ret\n",
    "\n",
    "def multinomial_basis(x, feature_num=10):\n",
    "    x = np.expand_dims(x, axis=1) # shape(N, 1)\n",
    "    feat = [x]\n",
    "    for i in range(2, feature_num+1):\n",
    "        feat.append(x**i)\n",
    "    ret = np.concatenate(feat, axis=1)\n",
    "    return ret\n",
    "\n",
    "def gaussian_basis(x, feature_num=10):\n",
    "    centers = np.linspace(0, 25, feature_num)\n",
    "    width = 1.0 * (centers[1] - centers[0])\n",
    "    x = np.expand_dims(x, axis=1)\n",
    "    x = np.concatenate([x]*feature_num, axis=1)\n",
    "    \n",
    "    out = (x-centers)/width\n",
    "    ret = np.exp(-0.5 * out ** 2)\n",
    "    return ret"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 返回一个训练好的模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(x_train, y_train):\n",
    "    \"\"\"训练模型，并返回从x到y的映射。\"\"\"\n",
    "    \n",
    "    basis_func = identity_basis\n",
    "    phi0 = np.expand_dims(np.ones_like(x_train), axis=1)\n",
    "    phi1 = basis_func(x_train)\n",
    "    phi = np.concatenate([phi0, phi1], axis=1)\n",
    "    \n",
    "    #==========\n",
    "    dim = phi.shape[1]\n",
    "    w = 0.05 * np.random.randn(dim)\n",
    "    for i in range(10000):\n",
    "        delta_w = np.dot(phi.T, y_train - np.dot(phi, w))\n",
    "        w += 0.001*(delta_w/np.mean(np.abs(delta_w)))\n",
    "    #==========\n",
    "    \n",
    "    def f(x):\n",
    "        phi0 = np.expand_dims(np.ones_like(x), axis=1)\n",
    "        phi1 = basis_func(x)\n",
    "        phi = np.concatenate([phi0, phi1], axis=1)\n",
    "        y = np.dot(phi, w)\n",
    "        return y\n",
    "        pass\n",
    "\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 评估结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(300,)\n",
      "(200,)\n",
      "训练集预测值与真实值的标准差：2.0\n",
      "预测值与真实值的标准差：2.2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xmc1XP7+PHXNUuFspQlS8l959amopFGSWTttpRufd2WCJ2ZkiitfiLSrYUkbTMqDNmjkrhLZgodKS1Kk7siJJQhItUs1++P9znmzDQzzXLOnJlzrufjcR4z53M+y/vT8LnOe7veoqoYY4yJXjHhLoAxxpjwskBgjDFRzgKBMcZEOQsExhgT5SwQGGNMlLNAYIwxUc4CgQk7ETlfRL4IdzkigYg0FJHfRSQ23GUx1YcFAlNpRGSbiFxceLuqfqCqZ4SjTIWJyEgRyfY9THeLyHIRSQx3uUpLVb9R1dqqmhvuspjqwwKBiVoiElfMR6+oam3gWCAdeK2Sr29MpbJAYMJORDqJyPaA99tEZJCIfCYiv4rIKyJSK+DzK0VkbcA39pYBnw0Tka0iskdENopIt4DPbhWRj0TkCRH5GRhZUrlUNQeYDZwsIseV8vpni8ga3/Vf85X9kcD7FJGhIvID8EwpzjdURL7zne8LEens295WRFaJyG8i8qOITPBtbyQi6g8yInKSiMwXkZ9FZIuI9A4490gReVVE0nzn/1xEEkr9hzMRwwKBqap6AJcDpwEtgVvBPWiBWUASUA9IAeaLSE3fcVuB84GjgIeAF0TkxIDzngt8CRwPjC6pACJSA+gJZAG/HOr6vv3fBJ4F6gIvAd0Knba+77NTAc8hzncG0A84R1XrAJcB23zneRJ4UlWPBP4OvFrMbbwEbAdOAv4F/McfTHyuBl4GjgbmA5NL+jcxkckCgamqJqnqDlX9GXgLaO3b3htIUdUVqpqrqs8B+4F2AKr6mu+4PFV9BdgMtA047w5VfUpVc1T1z2Ku3UNEdgN/+q73L1/t4FDXbwfE+cqerapvAJ8UOnce8KCq7vddv6Tz5QI1gWYiEq+q21R1q+882UBjETlWVX9X1Y8L34SINAA6AENVdZ+qrgVmADcH7Pahqi709Sk8D7Qq5t/ERDALBKaq+iHg971Abd/vpwL3+ppRdvse2A1w33gRkZ4BzSy7gRa4tn6/b0tx7VdV9WjgBGAD0Cbgs5KufxLwnRbM5Fj4ertUdV9pzqeqW4B7cE1YO0XkZRE5yXfc7cA/gE0islJEriziPk4CflbVPQHbvgZODnhf+N+5lvVdRB8LBKa6+RYYrapHB7wOV9WXRORU4Glcc0o938N8AyABx5c63a6q/oRrshkZ0LxU7PWB73H9CYHXa1D4tKW9H18ZXlTVDriAocBY3/bNqvpvXBPXWOB1ETmi0Ll3AHVFpE7AtobAd6X9NzDRwQKBqWzxIlIr4FXWb59PA8kicq44R4jIP30PuyNwD8tdACLSC1cjKDdV3QT8FxhSiut7cc05/UQkTkSuoWCzVJnuR0TOEJGLfP0f+3BNVbm+e7tJRI5T1Txgt+9cBYaMquq3wHLgUd+/dUtcTWJ2Rf5NTOSxQGAq20LcA83/GlmWg1V1Fa5dfTKuA3cLvo5kVd0IPI57IP8InAl8FIQyj8d17B5/iOsfAK7FPWx3AzcBC3Bt/mW+H1z/wBjgJ1wTzvHAfb7PLgc+F5HfcR3H1xdqcvL7N9AIVzt4E9c/sbiM928inNjCNMaEjoisAKar6jPhLosxxbEagTFBJCIXiEh9X9PQLbihr++Gu1zGlMRGBxgTXGfgxvTXxs1p+Jeqfh/eIhlTMmsaMsaYKGdNQ8YYE+WqRdPQscceq40aNQp3MYwxplr59NNPf1LV4w61X7UIBI0aNWLVqlXhLoYxxlQrIvJ1afYLWdOQiDQQkXQRyfRlNbzbt32kL5viWt+rS6jKYIwx5tBCWSPIAe5V1dW+WZefioh/IssTqvpYCK9tjDGmlEIWCHxD5r73/b5HRDIpmOzKGGNMFVApfQQi0gg4C1gBtMflYukJrMLVGn4p4hgP4AFo2LDhQefMzs5m+/bt7NtX1Kz6yFGrVi1OOeUU4uPjw10UY0yECvk8AhGpDSzFZVh8Q0ROwOVOUWAUcKKq3lbSORISErRwZ/FXX31FnTp1qFevHgWTPUYOVSUrK4s9e/Zw2mmnhbs4xphqRkQ+VdVDrjoX0nkEIhIPzAFm+xbpQFV/9C3AkYfLvHio7IxF2rdvX0QHAQARoV69ehFf6zHGhFcoRw0JMBPIVNUJAdsDlw3shssXX95rlL+A1UQ03KMxphheLzz6qPsZQqHsI2iPWxJvvYis9W27D/i3iLTGNQ1twy38YYwxJpDXC507w4EDUKMGLFkCiYkhuVTIagSq+qGqiqq2VNXWvtdCVb1ZVc/0bb+6uibk2r17N1OnTi3zcV26dGH37t2H3tEYE9Xy3n+fffv3Q26uCwYZGSG7luUaKqfiAkFubm4Re+dbuHAhRx99dKiKZYyJAJ9//jkdX32VgSIQG+tqBJ06hex61SLFRNB4vS6qdupU4SrWsGHD2Lp1K61btyY+Pp7atWtz4oknsnbtWjZu3EjXrl359ttv2bdvH3fffTcejwfIT5fx+++/c8UVV9ChQweWL1/OySefzLx58zjssMMqfp/GmGpp3759PPLII4wbN446derQe/hwOPzwoDyzSqSqVf7Vpk0bLWzjxo0HbSvR8uWqhx2mGhvrfi5fXrbjC/nqq6+0efPmqqqanp6uhx9+uH755Zd/fZ6VlaWqqnv37tXmzZvrTz/9pKqqp556qu7atUu/+uorjY2N1TVr1qiq6nXXXafPP/98kdcq870aY6qdJUuWaOPGjRXQm2++WXfu3FnhcwKrtBTP2OhpGsrIcO1sIWpva9u2bYGx/pMmTaJVq1a0a9eOb7/9ls2bNx90zGmnnUbr1q0BaNOmDdu2bQtqmYwxVd9PP/3ELbfcQufOnVFVFi9eTFpaGscdd8ikoUETPU1DnTq5djZ/D3yQ29uOOOKIv37PyMjgvffew+v1cvjhh9OpU6ci5wLUrFnzr99jY2P5888/g1omY6KG1wtpafDDD+79zz/DN9/AYYfBPfeAr2m2KlFV0tLSuPfee/n111+57777uP/++8PSPBw9gSAx0Q2/ClIfQZ06ddizZ0+Rn/36668cc8wxHH744WzatImPP/64QtcyxhTD64Vx42D+fMjLK3qfpCSYPRvGjAltO3sZbN68meTkZN5//33OO+88UlJSaNGiRdjKEz2BANx/BEH6D6FevXq0b9+eFi1acNhhh3HCCSf89dnll1/O9OnTadmyJWeccQbt2rULyjWNMQG8XrjwQti//9D7LlsG7dvD4MEwdmzoy1aMAwcOMG7cOB555BFq1arFtGnT8Hg8xMSEuZW+NB0J4X4FpbO4GoumezWm1P7zH1URVSjbq2vXCg8WKY8PPvhAmzVrpoBed911umPHjoN3Wr7c3VeQykcpO4ujq0ZgjKn+/MPA69Vz/X3+GkFMDHToAHXruj6C//0vv88g0Ny58M47kJ5eKU1Fv/zyC8OGDSM1NZWGDRuyYMEC/vnPfx68Y2oq3Hmna+KqWTOkM4kLs0BgjKk+Cj8sJ02CNWvcZz17Hvzg9Hph2DDXNBRo/37XuRzCB62q8uqrr3L33Xeza9cuBg4cyEMPPUTt2rUP3tnrhb593ahGf/kyMiwQGGNMAUU9LLOyYNq04o9JTISlS10AGTUKtm/P/2z1anfOEDxst23bRt++fXnnnXdo06YNCxcu5Oyzzy7+gLS0/PsCEAnpTOLComcegTGmeqvIw9LjgVdfdU1JfqtWuaRuQczsmZOTw2OPPUbz5s1ZtmwZEydOZMWKFSUHgaJcdVWljnCyQGCMqZ7K+rBMTHTNLZde6voT8vJcrWLkyKAEg5UrV5KQkMDgwYO5+OKLyczM5O677yY2NrbkA/3Xjo93wa1mTRgypMLlKQsLBMaY6qFnT/eNXsT9LM/DMjHRPfhr1swPBu+9V6GawZ49e+jfvz/nnnsuu3btYs6cOcydO5cGDRoc+mB/qumnn3blSUqqtE7sQBYIyqm8aagBJk6cyN69e4NcImMinP8b/ejRFetI9U8uvfji/GCwb59reiqjuXPn0rRpUyZPnkzfvn3ZuHEj1157bekXlMrIcLWS3FzIzoaGDcMy6c0CQTlZIDAmDBITYfjwij8s/TWD+Hj3XhVmzSp1rWD79u1069aNbt26UbduXZYvX87kyZM56qijylaO3bvzZ0Tn5bkhsWFgo4bKKTAN9SWXXMLxxx/Pq6++yv79++nWrRsPPfQQf/zxBz169GD79u3k5uYyYsQIfvzxR3bs2MGFF17IscceS3p6erhvxZjolJgIvXpBSooLBLm5h6xp5ObmMnXqVP7f//t/ZGdnM2bMGAYOHEi8P6CUhdcLjz+e/17EjYIKg4gIBPfccw9r16499I5l0Lp1ayZOnFjs52PGjGHDhg2sXbuWRYsW8frrr/PJJ5+gqlx99dUsW7aMXbt2cdJJJ/H2228DLgfRUUcdxYQJE0hPT+fYY48NapmNiShBXD+kWD17wnPPlSoZ5bp16/B4PHzyySdceumlTJs2jb/97W/lv3bhUVAxMZU6ZDRQRASCcFu0aBGLFi3irLPOAuD3339n8+bNnH/++QwaNIihQ4dy5ZVXcv7554e5pMZUE/48Qv4HdKg6UEuRjPKPP/7goYceYsKECdStW5fZs2fz73//u/T9AKVVyUNGA0VEICjpm3tlUFWGDx9OUlLSQZ99+umnLFy4kOHDh3PppZfywAMPhKGExlQzaWn5qSNCPQu4cDLK1FSYMwe6d+fdhg3p06cP27Zt4/bbb2fcuHHUrVs3ONft2ROeeSY/2FXykNFAEREIwiEwDfVll13GiBEjuPHGG6lduzbfffcd8fHx5OTkULduXW666SZq167Ns88+W+BYaxoyphhF5QiqDKmpkJTEj8A9ixbxMnDGGWewdOlSOnbsGNxrJSa6mk6om79KwQJBOQWmob7iiiu44YYbSPT9IWvXrs0LL7zAli1bGDx4MDExMcTHxzPNNxXe4/FwxRVXcOKJJ1pnsTGFeb0uKZxfXJz79lwJ8l5/nZnAEGAvMFKEYSkp1Ax2EPALYmr8CilNitJwvywNdfTcqzGanJyfXlrEva8EGzdu1A6+NYMvAM30p61u2za4aauDnGq6JNiaxcaYasfrdeP5Vd37GjVCXhvYt28fDzzwAK1ateLzrCxmXn456UAT/w6ffOKabiqahsLrhW7d4Pzz4f77g57nqCIsEBhjqo60NDfDFty4+l69Qtp0kp6eTqtWrRg1ahQ9evRg06ZN3PbOO0hycsEdDxwo18zjv3i9LpjMneuGjPrzHGVkVKT4QVOtA4H6vzVEsGi4R2MA97CcOTO/NhAfH7LaQFZWFr169eKiiy4iJyeHRYsW8cILL3D88ce7HXr2dPmIAqWmwtCh5btgWpoLJoHCOG+gsGobCGrVqkVWVlZEPyhVlaysLGrVqhXuohgTeoG1AYAuXYJeG1BVnn/+eZo0acILL7zAsGHDWL9+PZdccknBHf0jetq2zd+WlwfjxrmAUBapqQcfExsLU6ZUjY5iqvGooVNOOYXt27eza9eucBclpGrVqsUpp5wS7mIYU/nq1w/q6bZs2UJycjJLliyhXbt2pKamcuaZZxZ/QGIiTJzolr/05wMCGD8ezjyzdA/x1FTo06fg8c2awYwZVSYIANV31JAxJsKkpKjGxrqRQjVrBm1Uzf79+3X06NFaq1YtPfLII3XKlCmam5tb+hMMGZK/8L3/FR/vyluSlJT80U+Br0oaBaVa+lFDYX/Il+ZlgcCYCLd8uephh6nGxKjGxR36IVtKH330kTZv3lwB/de//qXfffdd+U6UkqLauHHBB3pcnAsSl15asLzLl7uHfUzMwUEgLq5Sho36lTYQVNumIWNMBElLc2sCqAYlC+fu3bsZPnw406dPp0GDBrz11ltceeWV5T+hx+Oagy64IL8fIyfH9RkALFrk0lLs3g1r1rjPCvdfVsUmIZ9q21lsjIkQqaluhS7/gzM2ttyjaVSVV199laZNm5KamsqAAQPYuHFjxYKAX2IiTJ7sZjrHxLhXoEWL3JyD7OyDg0B8fJUNAhDCzmIRaQCkAfWBPCBVVZ8UkbrAK0AjYBvQQ1V/CVU5jDFVmNcLd96Zn45ZBG67rVwPzK+//pq+ffuycOFCzj77bBYsWECbNm2CW15/zSAjw33799cIilKjhhv5VL++G45aRYMAhHbUUA5wr6quFpE6wKcishi4FViiqmNEZBgwDCjn4FxjTLWWkVFwRE058grl5OTw5JNP8sADDyAiTJgwgbvuuou4uBA93grnB3rssYL3ANC1q8smWoUf/oFCFghU9Xvge9/ve0QkEzgZuAbo5NvtOSADCwTGRKdOndzErf37XVPL5MlleniuWrUKj8fDmjVruPLKK5k8eTKnnnpq6Mpb2Nix7qHvryGsXQvdu7uaQzUiWgkTskSkEbAMaAF8o6pHB3z2i6oeU8QxHsAD0LBhwzZff/11yMtpjAmDcqxEtmfPHkaMGMFTTz3FCSecwKRJk+jevXvwF4up5kTkU1VNOOR+oQ4EIlIbWAqMVtU3RGR3aQJBoISEBF21alVIy2mMqR7mz59Pv3792L59O8nJyTz66KNlXzQ+SpQ2EIR01JCIxANzgNmq+oZv848icqLv8xOBnaEsgzGmivJ64dFHS52B87vvvqN79+5cc801HHXUUXz00UdMnTrVgkAQhHLUkAAzgUxVnRDw0XzgFmCM7+e8UJXBGFNFeb0uDbN/mcYlS4ptFsrNzWX69OkMHz6c7Oxs/vOf/zBo0CDi4+MrudCRK5Q1gvbAzcBFIrLW9+qCCwCXiMhm4BLfe2NMNMnIcEEgN9f9LCYd82effUb79u3p168f5557Lhs2bGD48OEWBIIslKOGPgSK67npHKrrGmOqgU6dXE3AXyMoNIFs7969PPzwwzz++OMcffTRPP/889x4443WGRwilmLCGBMet9zifhaabLVo0SKSk5P56quv6NWrF+PHj6devXphKmR0sEBgjKlcqanQr59rFqpZ868JZDt37mTAgAG8+OKL/OMf/yA9PZ1OVWThlkhngcAYU3n8KSVyctz7/fvR9HRmbdzI4MGD+eOPP3jwwQcZPnw4NQuvEGZCxgKBMabyFEopsSkmhqTXXmPZ2rV07NiRlJQUmjRpUvzxJiQs+6gxpvL4UkrsE+FBEVqqsv7rr5kxYwbp6ekWBMLEAoExpvIkJpIxfjyt6tXjYVWu+7//IzMzk9tvv52YwmmdTaWxf3ljTKXIysritttu48J+/ciuU4d3332X2bNnc8IJJ4S7aFHPAoExJqRUlRdeeIGmTZuSlpbG0KFD2bBhA5dddlm4i2Z8rLPYGBMyW7dupU+fPixevJhzzz2X9957j5YtW4a7WKYQCwTGmKDLzs7mscce4+GHHyY+Pp7J995L8jHHEPvHH+EumimCBQJjTFB5vV48Hg8bNmzg2muvZdItt3Dy9de7dBKjR5eYYM6Eh/URGGOC4tdff6Vv3760b9+e3bt3M2/ePObMmcPJn39eqgRzJnwsEBhjKkRVef3112natCkpKSn079+fjRs3cvXVV7sd6tVzy1DGxBSZYM6EnzUNGWPK7ZtvvuHOO+9kwYIFnHXWWcyfP5+EhIAFsbxeuOceVxuIiYGJE61ZqAqyQGBMVeNfw7dePVizBn74AX7+Gfbtg9tvrxILo+fk5PDUU08xYsQIVJXHH3+c/v37ExdX6JGSlubKrQoikJUVngKbElkgMKYqSU11Sdlyc93Ds7BPPoGtW2Hs2Movm8+nn36Kx+Nh9erVdOnShSlTptCoUaODd/R6YebM/PuIjbVmoSrK+giMqQq8XujTB/r2dZk5iwoCfuPHu31LudZvsPz+++8MGDCAtm3bsmPHDl555RUWLFhQdBAAVxvIzs5/36WLNQtVURYIjAk3//q9KSmuJnAoqm7fCy5wNYhKsGDBApo3b87EiRPp3bs3mZmZ9OjRo+QVw374oeD7+vVDW0hTbtY0ZEw4+TtT//wzf5sIxMXBgAHw22/ugbptG6xbl19TUHXftpOT4Z13YMiQkHzb3rFjB3fffTevv/46zZo148MPP6R9+/aHPjA1Fd56K/99fPxfC9CYKkhVq/yrTZs2akzEWb5cNT5e1T3W3Ss+XjU52X1W1P7JyaqxsQWPAbctJSVoRcvNzdWpU6fqkUceqTVr1tTRo0fr/v37y39fyclBK5spPWCVluIZazUCY8KlcBs6uFFB06YVvX9ionuddZbrIwhY4IXcXFc7gAqPKtqwYQMejwev10vnzp2ZPn06jRs3Lv0JMjIKNnHFxVltoIqzPgJjqorY2NI9MD0eFyxiYwtuV61QJ/Kff/7Jfffdx1lnncXmzZtJS0tj8eLFZQsC8NfiM8TEuCahKVOsk7iKs0BgTLj07Olm2oq4h/rUqaV/YHo88MEH0KxZwe15eTBuXJmLsnjxYlq0aMGjjz7KTTfdRGZmJjfffHPJncHFSUx0+YQeeQSWLq0S8x5MyURLGqZWRSQkJOiqVavCXQxjgsM/Ycw/pt7/e3m+NXu90KFDwWYiEZg+vVQP4J07dzJw4EBmz57N6aefTkpKChdeeGHZy2GqJBH5VFUTDrWf9REYU5n8Q0UPHHC1gSVLYPjw8p8vMdE1EwX2GfibiKDYYKCqPPPMMwwePJg9e/YwYsQI7rvvPmrVqlX+skDBIGfNQdWGBQJjKlNgygV/Js6KPjD9D/vAYJCX5yannXnmQef/4osvSEpKYunSpXTo0IGUlBSaFW5iKo+igpwFg2rB+giMqSxeLzzzTGhSLvg7kAPb9HNzC/QX7N+/n4ceeoiWLVuybt06UlNTWbp0aXCCALigZummqyULBMZUlowMlz4C3AP7ttuC+43Z44Frrim4bd48GDqUZcuW0bp1a0aOHMm1115LZmYmvXv3JiYmiI+ATp1cTSA21tJNVzPWWWxMZamMphOvF84//69x/D8DQ4CZQKNGjZg2bRqXX355cK/pv64/Y2pWlvURVBHWWWxMVeJ/UE6cGNoHZWIiTJ2KJiXxEjAAyAIGAw8OHMgRoQoCF16YH+DS0y0IVDMWCIwJNX9q6bw8N9EqxJ2oX158MX0aNWLRtm2cA/wXaA0ud1FCQvCvPW4c7N/vft+/33WIWyCoVkLWRyAis0Rkp4hsCNg2UkS+E5G1vleXUF3fmCrB64V+/VzfQF6ee1CGqBM1OzubsWPH0qJFC5b/9BOTOnfGiy8IwEGdx0GRmur6IUy1FsrO4meBouqhT6hqa99rYQivb0z4Fc67ExMTkk7UFStWkJCQwLBhw7j88svJzMzkrvfeI7Zr14I7zpsXvNTV/iAX2M9Y2jQZpkoJWSBQ1WW4vipjole9eu7hGKK8O7/99hv9+vUjMTGRrKws5s6dyxtvvMEpp5zidhgypGBOIlX38A7GojZpafmjoMDdY1nSZJiqozQpSsv7AhoBGwLejwS2AZ8Bs4BjSjjWA6wCVjVs2DDo6VmNCbnly1Vr1gxJmui8vDydM2eOnnTSSSoi2r9/f/3tt9+K3jklRTUmpmBa6LZti051XVpDhqiKhCwNtgkOSpmGurLnEUwD/o5rtvweeLy4HVU1VVUTVDXhuOOOq6zyGRM8aWn5nai5uW4h+iD49ttv6dq1K927d+e4445jxYoVPPnkk9SpU6foA4rKVvrJJ26YaXmaiYYOdX0NgU1CvXtbcrlqrFIDgar+qKq5qpoHPA20rczrG1Od5ebmMnHiRJo2bcp7773H+PHjWbVqFeecc86hD/Z43MO64AldWoqyBAOvFx57rOC2mBjrF6jmKjUQiMiJAW+7ARuK29eYai8wzXSNGhV6WK5Zs4Zzzz2XAQMG0LFjRz7//HMGDRpEXFwZRoD37OmGrwby5yQqTZ+Bf1nNwEynAIMGWb9ANRfK4aMvAV7gDBHZLiK3A+NEZL2IfAZciJvvYkzk8U8ge+opGD263Mnlfv/9d+69914SEhLYvn07L7/8Mm+//TaNGjUqe5kSE91kr8IjiUozrHToUJfu+pNP8rfFxLjO6LFjy14WU6VYigljgi1IqSTefvtt+vbtyzfffIPH42HMmDEcc8wxwSljt24wd27+e//COIGznv3B7PPPYfbsg8+RnFz8spqmSrAUE8aES1FZOMsQCL7//nvuvvtuXnvtNZo1a8YHH3xAhw4dglvGIUPgrbfy5zjk5blhpXl57pv+6afD5s3u88JNQVDhpi5TtVj2UWOCrV491y8QE1OmLJx5eXlMnz6dpk2bMn/+fEaNGsWaNWuCHwTgr5xExMW5csbGuod+bi5kZ8PGje5nUUGgY8fgrKNgqgyrERgTTIEdqrGxLslcKR6YGzZsICkpieXLl3PhhRcyffp0/vGPf4S2rB6PW7jGnzX0nnvgzz+L3lfEvQYNsj6BCGSBwJhgClyBTMS1uZfgzz//5JFHHmHcuHEcddRRPPvss/Ts2bN8i8aXR2JifqA680zXaRzYd9C1K1xxhaWWjnAWCIwJFq8XZs3Kn2gVF1dis9CSJUtITk5my5Yt3HLLLTz22GMce+yxlVPWoiQmwptvunkFc+ZA9+42SSxKWCAwJljS0ly7OrjaQK9eRX6D3rVrF/feey/PP/88jRs3ZsmSJVx00UWVXNgSeDwWAKKMdRYbEwxeL8ycmV8biI8/aFSNqvLss8/StGlTXn75Ze6//37Wr19ftYKAiUpWIzAmGAJrAwBduhSoDfzvf/8jOTmZ9PR02rdvT0pKCs2bNw9DQY05mNUIjAmGH34o+L5+fQD279/PqFGjaNmyJatXr2b69OksW7bMgoCpUqxGYExFpaa6yVl+cXHQsycffPABSUlJZGZm0qNHDyZOnMiJJ55Y/HmMCROrERhTEV6vW4/YP0NXhF9uvhnPM8/QsWNH9u7dy9tvv83YjQEvAAAcPUlEQVQrr7xiQcBUWVYjMKYiApaiVOCVmBjumTePn379lUGDBjFy5EiOOOKIsBbRmEOxQGBMRdSrB6p8BfQF3s3NJeFvf+Pdp5+mdevWhzramCrBAoExFZC9ahUTgQeBWODJK6/kzrlziQ1cDcyYKs4CgTHl9MnMmXiefpp1wDXAUzVq0OC++wouCWlMNWCdxcaU0W+//cZdd91FuzvuYBfwBjAXaFBo7oAx1cUhA4GI9BORIK2GYUz19uabb9KsWTOmTJnCnaedRiZuzVXgr7kDxlQ3pakR1AdWisirInK5VFpaRGOqju3bt9O1a1euvfZa6tWrhzc1lad27OBI/w5FpJQwpro4ZCBQ1fuB04GZwK3AZhH5j4j8PcRlMybscnNzmTRpEk2bNmXRokWMHTuWVatWce6uXZCT43YSgdtvt2YhU22VqrNYVVVEfgB+AHKAY4DXRWSxqg4JZQGNCZe1a9fi8XhYuXIll112GdOmTeO0005zH3bq5FYf869LbLUBU40dMhCISH/gFuAnYAYwWFWzRSQG2AxYIDAR5Y8//mDkyJE88cQT1KtXjxdffJHrr7++4GIxiYluUfqMDFuwxVR7pakRHAtcq6pfB25U1TwRuTI0xTImPN555x369u3Ltm3buOOOOxg7dix169YteufA1b2MqcZK00fwQOEgEPBZZvCLZEzl++GHH7j++uvp0qULtWrVYtmyZTz99NNFBwGvFx591P00JgLYhDIT1fLy8pgxYwZDhw5l7969PPzwwwwZMoSaNWsWfYDXC5075/cNLFlitQJT7dmEMhO1Nm7cSMeOHUlKSqJ169asX7+eESNGFB8EIH9x+txcFwwyMiqtvMaEitUITOilprplHGvVgsCmlvr13WibSv5GvW/fPkaPHs3YsWOpU6cOzzzzDLfccguHnCLj9cIzz+QvRxkbW+Li9MZUFxYITGh4ve7b88cfw9q1xe83YwbccUelBYT333+f5ORkNm/ezM0338zjjz/OcccdV7qDMzIKzh247TZrFjIRwZqGTPB5vXDBBTB9eslBANyDNSXF7Z+aGrIi/fTTT9x666107tyZvLw8Fi9eTFpaWumDALiU0zEx7lWrls0dMBHDAoEJLq/XzbINXMj9UFTd/snJ0K1bUEfjqCrPPfccTZo0Yfbs2dx3332sX7+eiy++uGwnSk11K5Hl5LhAMHGi1QZMxLCmIRM8Xq9rMz9w4ODPRKBVK2jUqOD2t97KX+ZRFebOhbffhqVLK/yg3bx5M8nJybz//vucd955pKSk0KJFi7KfyOuFfv3ym4Xy8iArq0JlM6YqsUBggicj4+CagAhccw0MGVL0gz01Ffr0cQ9Xv+xsGDbMBYNyOHDgAOPGjeORRx6hVq1aTJs2DY/HQ0xMOSvAgX0D4GoE1klsIog1DZng8Hrhm28KLsoSG+v6Cd58s/hv9x4PTJvmAkagZcvK1Wfw0UcfcdZZZzFixAiuvvpqMjMzSU5OLn8QANi9O3+kEMDAgdYsZCJKyGoEIjILuBLYqaotfNvqAq8AjYBtQA9V/SVUZTCVJHCSVVwcXHll2YaGejzuZ1JSwe3jx8OZZ5bqHLt372bo0KGkpqbSsGFDFixYwD//+c9y3EwhXi9MmJD/PiYGjj664uc1pgoJZY3gWeDyQtuGAUtU9XRgie+9qe4CJ1nl5EDbtu5bflm+NXs8rvko0JYt0LFjiTUDVeWVV16hSZMmzJgxg4EDB/L5558HJwiAaxYKbLayuQMmAoUsEKjqMuDnQpuvAZ7z/f4c0DVU1zeVJDUVnn46OJOsxo51Q0kbN87flpMDffsWOZJo27Zt/POf/+T666/nlFNOYeXKlTz++OPUrl27fNcvin/IqIir7UyebM1CJvKoasheuCagDQHvdxf6/JcSjvUAq4BVDRs2VFMFLV+uGhen6sKAqohqcnJwzhsbm39eKHDe7OxsHT9+vB5++OF6xBFH6BNPPKHZ2dkVv25R5TjsMNWYGNX4eNWUlOBfw5gQAlZpKZ7VVbazWFVTVTVBVRPKNOnHVJ7CzSZxccGZZJWYCFddVXDb6tXg9bJy5UrOOeccBg8ezMUXX0xmZib33HMPcXEh6O7yN3nl5dmQURPRKjsQ/CgiJwL4fu6s5OubYPGPEoqLc00nwW42GTLEZff02bNyJXeffz7t2rVj586dzJkzh7lz59KgQYPgXK8wyytkokhlB4L5uNXO8P2cV8nXN8HgHyX09NOu7dzjccM9/aN/giEx0dU4Lr2UeSI0U+Wp3Fz6tG3Lxo0bufbaaw+dJK4i0tLyJ8ZZXiET4UIWCETkJcALnCEi20XkdmAMcImIbAYu8b031U1GhntI+kcJNWwYkofkdw0bcu3+/XRV5RhgeXw8k1u35qiNG4N+rQK8Xpg1K782YGsSmwgXylFD/1bVE1U1XlVPUdWZqpqlqp1V9XTfz8Kjikx1EJh8rUaNoDeZ5ObmMnnyZJo2bco7K1Ywpk8fPu3dm3YiblRRhw4wdGhQr1lAWlr+DGkR6NXLagMmolXZzmJTRXm9cNdd7kEpEvTka+vWreO8887jrrvuol27dmzYsIGhU6cSf9pp7pqqruN23LjQZCstPBw2Pt5qAybiWSAwZRPYdp6bC2vWBOW0e/fuZejQobRp04avvvqK2bNn89///pe///3vbodOnQ5OQzFzZlCu/Rd/3iN/EjywvgETFSwQmLL54Yegn/K///0vLVq0YNy4cdx6661s2rSJG264oWBncGIiDBpU8MCVK4OXttrrdWmmQzEc1pgqzgKBKb3UVJc22q+CzSY//vgjN9xwA5dffjk1atRg6dKlzJgxg7qBy1kGGjsWugZMRvenre7YseLBIC2tYE0gJgamTLHagIkKFghM6fi/MfsfliJuAZpyPCjz8vKYMWMGTZo0Yc6cOYwcOZJ169bRsWPHQx88ZIj7ph4oJ8f1GZSX1+uamQLnDEybFtzhsMZUYRYITOkEaRZxZmYmnTp1onfv3rRq1Yp169bx4IMPUrNmzdKdIDHRfVMvnFZ63rzydx6PG1dwHYWrrrIgYKKKBQJTOvXquW/K5Uy+tm/fPh544AFatWrFhg0bmDlzJunp6TRp0qTsZSlqDQNV19Fb1mDg9RZs7gKXQtuYKGKBwBxa4JBRf9t5Gb4xp6en06pVK0aNGkWPHj3YtGkTt912W8VmBns8btGbwJpBXp4LBmXpQM7IKLjoTGysdRCbqGOBwBxaOYeMZmVl0atXLy666CJycnJYtGgRL7zwAscff3xwyuWvGRQOBnPnuklnxdUOvF549NH8NZZr1szPlzR1qnUQm6hjaxaboFNVXnjhBQYOHMju3bsZNmwYI0aM4PDDDw/+xfw1k8LrHufluXUMtm6FtWuhe3e32llamksml5PjZkUvWeJeGRkuKFgQMFHIAoEpmb+JJT4+/+FZQtPJli1b6NOnD++99x7t2rUjNTWVM888M7Rl9AeD5OSCzTy5ufmjiRYtyr8H/z4HDrgAMHy4BQAT1SwQmOIVXos4KanYdYgPHDjAY489xqhRo6hRowZTpkyp+KLxZVFUzUCkYGAIHBkkEpI8ScZURxYITPEyMmD/fvdgVS02y+jy5ctJSkpiw4YNdO/enUmTJnHSSSdVfnk9nvzmH4Ajjyw4vyA+3t1LbKxLHVFMUDMm2lggMMXbvTv/23VenhtCWuDj3QwfPpzp06fToEED5s+fz1WFVxarbImJBR/uf/87zJmT30dgfQHGHMQCgSma1wtPPJH/XuSvpRpVlddff53+/fuzc+dOBgwYwMMPPxzcReODxeMpONTVAoAxB7FAYIqWluY6Vv18SzV+/fXX3Hnnnbz99tucffbZLFiwgDZt2oSvnMaYCrN5BOZghVfoio0lZ9IkJni9NGvWjIyMDCZMmMCKFSssCBgTAaxGYA6WkVEgudynXbvimTGD1atXc+WVVzJ58mROPfXUsBbRGBM8FgjMwTp1gho12LN/Pw+IMOnNNznhhBN47bXX6N69e2gXjTfGVDprGjIHS0xk/siRNK9Thyfz8khKSiIzM5N//etfFgSMiUBWIzAFfPfdd/Tv35833niDFi1a8EpqKok20saYiGY1AgNAbm4uU6ZMoWnTpixcuJD//Oc/rF692oKAMVHAagSGzz77DI/Hw4oVK7j44ouZPn16/qLxxpiIZzWCKLZ3716GDRtGmzZt2Lp1K88//zyLFi2yIGBMlLEaQZRatGgRffr04csvv6RXr16MHz+eeoVSSBhjooPVCKLMzp07ufHGG7nsssuIi4sjPT2dWbNmuSAQuGCLMSZqWI0gSqgqs2bNYvDgwfzxxx88+OCDDBs2jFq1arkdUlPhzjtdcrmaNd1iLdZRbExUsBpBFNi0aROdOnXijjvu4Mwzz2TdunWMHDkyPwh4vdCvn8stlJfnUk9nZIS1zMaYymOBIILt37+fkSNH0qpVK9avX8+MGTNIT0+nSZMmBXfMyCiYYC4mxhZsMSaKWNNQhFq6dClJSUl88cUX3HDDDUyYMIETTjih6J3r1Su4ktfAgdYsZEwUsRpBhPn555+5/fbb6dSpEwcOHODdd99l9uzZxQcBgDVr8n+PiYGjjw59QY0xVUZYagQisg3YA+QCOaqaEI5yRBJV5cUXX2TAgAH8/PPPDB06lAceeIDDDz+85ANTU+Hpp/Pfx8dbs5AxUSacTUMXqupPYbx+xNi6dSt9+vRh8eLFnHvuubz33nu0bNny0Af6O4n9KacBevWyZiFjoow1DVVj2dnZjBkzhhYtWvDxxx8zefJkPvroo9IFATi4kzguzi3oboyJKuGqESiwSEQUSFHV1MI7iIgH8AA0bNiwkotX9X388cd4PB7Wr1/Ptddey6RJkzj55JPLdpLdu62T2BgTtkDQXlV3iMjxwGIR2aSqywJ38AWHVICEhAQt6iRh4fW69Xx/+MG9r1/ffYuupAfor7/+yn333ce0adM4+eSTmTdvHldffXXZT1TU4vTWSWxMVApLIFDVHb6fO0XkTaAtsKzko8LMHwBmzoTs7IKfzZqVPwErI8N1tgY5MKgqb7zxBnfddRc//vgj/fv3Z9SoUdSpU6d8JyzcLORbnN4YE30qPRCIyBFAjKru8f1+KfBwZZejTLxe6NwZ9u0r2JTil53tgsRzz7l9RGDQIBg7NiiX/+abb+jXrx9vvfUWrVu3Zv78+SQkVHCglTULGWN8wtFZfALwoYisAz4B3lbVd8NQjtJJTXVNP8UFAXBDLiF/n7w8GDcOLrigQgnccnJymDhxIs2aNWPJkiU89thjrFy5suJBwOuFCRPy31uzkDFRrdJrBKr6JdCqsq9bLqmpkJSU/z4mxj30r7gif5u/j8C/f2CwWLYMOnSAadPA4ynTpVevXk3v3r1ZvXo1Xbp0YcqUKTRq1Kj89xIoLe3g0ULWLGRM9FLVKv9q06aNVrrly1UbN1Z1j3b3atzYbS/OkCEF9/e/YmJKPi7Anj17dMCAARoTE6P169fXV155RfPy8oJ0U+rKEReXX7bYWNWUlOCd3xhTZQCrtBTPWJtHUBR/n8DWrQW3Dx5ccjv62LEwZIhragmUlwd33HHIZqIFCxbQvHlznnjiCXr37k1mZiY9evRACp+vIgrXBq66qsy1FWNMZLFAUJS0tPz2/pgYaNwYUlJK98AcOxamTz84GGzcWGyfwffff891113HVVddRe3atfnwww+ZPn06R1dGu339+qG/hjGmSrNAUJg/946/rT8uzgWGsnxr9niKDgbZ2a4T2ScvL49p06bRpEkT3nrrLR555BHWrFlD+/btg3AjxejZ0y08I+J+2kxiY6KepaEOlJoKffq4phxwD8vbbivfsEp/4EhOLtiBPG8epKay4bzz8Hg8eL1eOnfuzPTp02ncuHHF76EkXq+bPzBpEmRlhWS+gzGm+rFA4OdPwOYPAlDx3DtFBIM/VRmVnMz4mBiOPuYY0tLSuOmmm4LbD1AUf7/HgQNQo4YtRWmM+YsFAr+MjIJZOGNiYPLkij8s/cGgTx/ey8sjGdiqyq2NGzP+ww859thjK3b+0grs9zhwwN2vBQJjDNZHkK9ePZdmwT9XoBxj/4uzq1s3bj7nHC7B/YO/DzyzZQvHvvFGUM5/SF6vS43hb6KydBLGmAAWCMA9KO+5x9UI/DWBIAQBVeWZZ56hSZMmvLJ6NSPatOEz4EJw1+rTx/VLhNq4cQXzI3XpYrUBY8xfLBBAfrNJXp771pyVVeFTfvHFF1x00UXcdtttNGvWjLVr1/LwU09Ry5+OAtz1+vWrUBqKQ/J64a23Cm6zIaPGmAAWCAoPF61gs8n+/ft5+OGHadmyJWvXriU1NZWlS5fSrFkz9y188mRX6/DLzc3PXBoKGRkFRy3FxtqQUWNMAdHdWez1wp135ncSV2S4KLBs2TKSkpLYtGkT119/PU888QT1C3/79jc5+ZeIrFkzdO31Xi98843r88jOdkEgGB3gxpiIEt01gsIjhco5XPTnn3/mjjvu4IILLmDfvn288847vPTSSwcHAT+PB5YuhUceccM416+Hyy4Lbn+B1+sCTEqKu0f/NS2dhDGmkOiuERTOyT9gQJm+LasqL730EgMGDCArK4vBgwfz4IMPcsQRRxz64MRE9wrMcLpokctvFIx1DNLS3DBRyM8tZDUBY0wRordGUMGlGr/88kuuuOIKbrzxRk499VRWrVrFuHHjShcEAs2ZU/D9+PGh7Tw2xphCojcQFG4WKmUncXZ2NmPHjqVFixZ89NFHTJo0Ca/XS+vWrctXju7dC75XdUNZKxIM/MfGx1tOIWPMIUVvICg8gWzKlEM2naxYsYKEhASGDRvGZZddRmZmJnfddRexsbHlL4fHAx07Ftz2yScuKJUnGKSmuiynqakuCCQlQXq6NQsZY4oVnYEgNdWNFsrJKdUEst9++41+/fqRmJhIVlYWb775Jm+++SannHJKcMozZoz71h7owIGy1wz8SfOys90chZwcaNjQgoAxpkTRFwj8yeVycvLXFy5mApmq8sYbb9C0aVOmTp3KXXfdRWZmJl27dg1umRIT3bf2tm0Lbv/kk/xv94cydKhLbheYNM9SSRhjSiH6AkFRyeWKeFh+++23dO3ale7du3PcccexYsUKnnzySerUqROaciUmwsSJB9cMsrPdA75bt+JrBzfd5NJIBI6AClbSPGNMxIu+4aP+vgEocoJVbm4ukydP5v777ycvL4/x48dzzz33EBdXCf9U/ppBWpqb7ewPWKowd65by+CGG+CPP2DHDjj9dFi9GjIzC55HJKhJ84wxEa40CxuH+xW0xeuXL1c97DC3mHxc3EGLtq9evVoTEhIU0CuuuEK/+uqr4Fy3PFJSXDn9i8yX5TVkSPjKbYypMrDF64tQTHK5P/74g0GDBnHOOefw7bff8vLLL/P222/TqFGj8JXV43Hf6mPK+Ce68cbgTEgzxkSN6AkExeTkX7hwIc2bN+fxxx/n9ttvJzMzk//7v/8L/YphpeEPBqUpS7NmLp3ECy+EvlzGmIgSPX0EaWkFcvJ/f9FF3P3EE7z22ms0a9aMDz74gA4dOoSxgMXweODMM11n8I4drmP7f//L7yPYtctNSrP+AGNMOUVPIPjhBwDygFRg2Pvvs0+EUaNGMWTIEGrUqBHW4pUoMRHefDPcpTDGRKjoCAReLyxcyOeAB1gOXNiiBdNfeol//OMfYS6cMcaEV1QEgj8XL2Z0djbjgCOBZy+8kJ5LllSNfgBjjAmziA8EH330EbdOmcIWVXoCj9eqxbGjR5euA9YYY6JAxAcCXb+emF27eE+EznFx8OSTNtvWGGMCRHwg6PDLL3wuQlxeXol5hYwxJlpFfCCgXj3i/JOyatSwJGzGGFNIWCaUicjlIvKFiGwRkWEhu5DX61I55+W5CWQTJ1qzkDHGFFLpgUBEYoEpwBVAM+DfItIsJBcLTClhzULGGFOkcNQI2gJbVPVLVT0AvAxcE/SreL0wa1Z+Som4OGsWMsaYIoQjEJwMfBvwfrtvWwEi4hGRVSKyateuXWW/SuC6AyLQq5c1CxljTBHCEQiKGsCvB21QTVXVBFVNOO6448p+lU6dXOdwbCzUqmWLtxtjTDHCMWpoO9Ag4P0pwI6gXyUxEZYscTWDTp2sNmCMMcUIRyBYCZwuIqcB3wHXAzeE5EqJiRYAjDHmECo9EKhqjoj0A/4LxAKzVPXzyi6HMcYYJywTylR1IbAwHNc2xhhTUPSsUGaMMaZIFgiMMSbKWSAwxpgoZ4HAGGOinKgeNJeryhGRXcDX5Tz8WOCnIBanOrB7jg52z9GhIvd8qqoeckZutQgEFSEiq1Q1IdzlqEx2z9HB7jk6VMY9W9OQMcZEOQsExhgT5aIhEKSGuwBhYPccHeyeo0PI7zni+wiMMcaULBpqBMYYY0pggcAYY6JcRAcCEblcRL4QkS0iMizc5akMIrJNRNaLyFoRWRXu8oSCiMwSkZ0isiFgW10RWSwim30/jwlnGYOtmHseKSLf+f7Wa0WkSzjLGGwi0kBE0kUkU0Q+F5G7fdsj9m9dwj2H9G8dsX0EIhIL/A+4BLcYzkrg36q6MawFCzER2QYkqGrETroRkY7A70CaqrbwbRsH/KyqY3xB/xhVHRrOcgZTMfc8EvhdVR8LZ9lCRUROBE5U1dUiUgf4FOgK3EqE/q1LuOcehPBvHck1grbAFlX9UlUPAC8D14S5TCYIVHUZ8HOhzdcAz/l+fw73P0/EKOaeI5qqfq+qq32/7wEyceubR+zfuoR7DqlIDgQnA98GvN9OJfyDVgEKLBKRT0XEE+7CVKITVPV7cP8zAceHuTyVpZ+IfOZrOoqYJpLCRKQRcBawgij5Wxe6Zwjh3zqSA4EUsS0y28EKaq+qZwNXAHf6mhRMZJoG/B1oDXwPPB7e4oSGiNQG5gD3qOpv4S5PZSjinkP6t47kQLAdaBDw/hRgR5jKUmlUdYfv507gTVwTWTT40de+6m9n3Rnm8oScqv6oqrmqmgc8TQT+rUUkHvdAnK2qb/g2R/Tfuqh7DvXfOpIDwUrgdBE5TURqANcD88NcppASkSN8HUyIyBHApcCGko+KGPOBW3y/3wLMC2NZKoX/YejTjQj7W4uIADOBTFWdEPBRxP6ti7vnUP+tI3bUEIBviNVEIBaYpaqjw1ykkBKRv+FqAeDWo34xEu9ZRF4COuHS8/4IPAjMBV4FGgLfANepasR0rhZzz51wTQUKbAOS/G3nkUBEOgAfAOuBPN/m+3Bt5hH5ty7hnv9NCP/WER0IjDHGHFokNw0ZY4wpBQsExhgT5SwQGGNMlLNAYIwxUc4CgTHGRDkLBMYYE+UsEBhjTJSzQGBMOYjIOb4EYLV8M7o/F5EW4S6XMeVhE8qMKScReQSoBRwGbFfVR8NcJGPKxQKBMeXky2G1EtgHnKequWEukjHlYk1DxpRfXaA2UAdXMzCmWrIagTHlJCLzcSvfnYZbXrBfmItkTLnEhbsAxlRHItITyFHVF33rYy8XkYtU9f1wl82YsrIagTHGRDnrIzDGmChngcAYY6KcBQJjjIlyFgiMMSbKWSAwxpgoZ4HAGGOinAUCY4yJcv8fROvYxwbT2GQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def evaluate(ys, ys_pred):\n",
    "    \"\"\"评估模型。\"\"\"\n",
    "    std = np.sqrt(np.mean(np.abs(ys - ys_pred) ** 2))\n",
    "    return std\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    train_file = 'train.txt'\n",
    "    test_file = 'test.txt'\n",
    "    x_train, y_train = load_data(train_file)\n",
    "    x_test, y_test = load_data(test_file)\n",
    "    print(x_train.shape)\n",
    "    print(x_test.shape)\n",
    "\n",
    "    f = main(x_train, y_train)\n",
    "\n",
    "    y_train_pred = f(x_train)\n",
    "    std = evaluate(y_train, y_train_pred)\n",
    "    print('训练集预测值与真实值的标准差：{:.1f}'.format(std))\n",
    "    \n",
    "    y_test_pred = f(x_test)\n",
    "    std = evaluate(y_test, y_test_pred)\n",
    "    print('预测值与真实值的标准差：{:.1f}'.format(std))\n",
    "\n",
    "    plt.plot(x_train, y_train, 'ro', markersize=3)\n",
    "    plt.plot(x_test, y_test_pred, 'k')\n",
    "    plt.xlabel('x')\n",
    "    plt.ylabel('y')\n",
    "    plt.title('Linear Regression')\n",
    "    plt.legend(['train', 'test', 'pred'])\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
